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1. Introduction 

Oh 

Parker (1966) first identified the instability of a gravitationally stratified gaseous disk 
\Q in a magneto-hydrostatic equilibrium, like the Galaxy, in response to perturbations due to 

i — i magnetic buoyancy, for perturbations that occur in the magnetic field lines that lie par- 

allel to the disk plane. This phenomenon is referred to as magnetic buoyancy instability 
(Hughes & Proctor 1988; Tajima & Shibata 1997) or Parker instability (Parker 1966; 1979). 
Solar magnetic activities such as sunspots, which is caused by the emergence of magnetic 
flux tubes from the interior of the sun into the solar atmosphere (Zwaan 1985, 1987), reflect 
this instability. The mushroom shape of the hydrogen cloud GW 123.4-1.5 (Baek, Ku- 
doh & Tomisaka, 2008) also suggested magnetic flotation and, hence, Parker instability. In 
addition, giant dense CO molecular loops close to the Galactic center has been ascribed to 
the instability (Fukui et al., 2006). As is expected, gas aggregating at the foot points of 
(-\| the rising magnetic loops eventually collapses, becoming a giant cloud and a region that 

hosts stellar formations. On the other hand, for models galactic dynamos, Hanasz et al. 
(2004) found that if the effect of cosmic rays is incorporated, the efficiency of sustaining the 
galactic magnetic field can be improved. Otmianowska-Mazur et al. (2009) suggested that 
the extended, X-shaped magnetic halo structures observed in some edge-on galaxies could 
be attributed to the cosmic-ray driven dynamo. Actually, based on his instability, Parker 
(1992) introduced the idea of cosmic-ray driven dynamo to maintain the Galactic magnetic 
X field, because the cosmic-ray pressure functions just like the magnetic pressure, capable of 

overpowering the gas pressure inside the magnetic flux tube and making it easier for the gas 
to rise. 
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Cosmic rays are a major component in the interstellar medium (ISM) in galaxies, whose 
energy density is comparable to the kinetic and magnetic energy densities of thermal plasma 
gas (Ferriere 2001). As high energy particles endowed with an energy of up to 10 14-15 eV, 
of which, 90% are protons, as long as their gyroradius is significantly smaller than the 
characteristic spatial scales of the magnetic field, cosmic-ray particles only propagate along 
the magnetic field lines. When the gyroradius is comparable to the scale of magnetic field 
variation, cosmic rays interact strongly with the field through gyroresonant scattering by the 
magnetic irregularities. Although the velocity of cosmic-ray particles is close to the speed of 
light, the bulk motion of cosmic rays is diffusive and the bulk speed is of the order of Alfven 
speed (Cesarsky 1980). 

The morphology of Parker instability exhibits two modes, i.e. the undular mode and the 
interchange mode. The undular mode, also called Parker instability, is excited by perturba- 
tions along the magnetic field lines (the wavenumber vector k\ \ of the perturbation parallel 
to the magnetic field B), where the falling gas creates a magnetic buoyancy greater than the 
restoring magnetic tension. The interchange mode, known as flute instability or magnetic 
Ray leigh- Taylor instability (Kruskal & Schwarzchild 1954), occurs for shorter- wavelength 
perturbations with fc_L perpendicular to B, capable of causing two straight flux tubes to 
interchange and ultimately reducing the potential energy in the system. The linear growth 
rate of the interchange mode generally exceeds that of the undular mode owing to a rapid 
growth of short wavelengths. However, in the nonlinear stage, the undular mode often dom- 
inates (Matsumoto et al, 1993; Tajima & Shibata 1997). Thus, the undular mode (Parker 
instability) is more important than the interchange mode in astrophysical problems. 

Baierlein (1983) and Matsumoto et al. (1988) performed the first one- and two-dimensional 
(2D), pure MHD simulations of Parker instability, respectively. In the nonlinear stage, the 
gas condensates to form giant clouds; in addition, a shock wave appears in the flow along the 
rising magnetic loop. By applying the simulation results of Matsumoto et al. (1988) to the 
solar atmosphere, Shibata et al. (1989b; 1990a) demonstrated that the emerging magnetic 
loop still expands self-similarly during the nonlinar evolution in two dimensions. Kamaya et 
al. (1996) adopted the supernova explosion as a perturbation in the ISM to trigger nonlinear 
instability. Earlier, Nozawa (1992) examined the instability deeper in the convectively un- 
stable layer of the solar atmosphere; when considering how magnetic shear affects the flow 
(Hanawa et al. 1992, Nozawa 2005), although the interchange mode is stabilized, a large 
unstable thin structure may still arise. Shantanu et al. (1997) and Kim et al. (2000) made 
a further application of Parker instability to the Galactic disk without cosmic rays. 

Parker's original analysis showed that the instability has a maximum growth rate for 
non-zero /c_L, although non-zero k\ | is the main cause for the instability. In the 3D simulations 
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for both solar and galactic problems given in Matsumoto & Shibata (1992) and Matsumoto 
et al. (1993), previous 2D results of cloud formation, presence of shock wave, and self-similar 
evolution are confirmed. However, the spatial and temporal scale in these studies depends on 
the k-L. If a larger fc_L is applied, the magnetic loop tends to have a thinner structure and a 
horizontal expansion, which would suppress the upward expansion (Parker 1979). Similarly 
thin slice structures have been found in other 3D simulations of Kim et al. (1998, 2001, 
2001) and Hanasz et al. (2002). 

This study describes the evolution of Parker instability undergoing cosmic-ray diffusion 
using 2D simulations. Effects of the adiabatic index and the initial hydrostatic condition 
other than an isothermal temperature and a uniform density profile are examined. Two 
initial profiles are used, i.e., the hyperbolic tangent temperature model, which is often used 
in the solar atmosphere, and the exponential density model, which is appropriate for Galactic 
problems. Various perturbation results are explored. The physical parameters and initial 
conditions invoked are appropriate for a galactic ISM in the solar neighborhood. Section 2-6 
details the governing equations, numerical algorithms, normalization, initial conditions, and 
grid setup. Section 7 presents the results. Section 8 draws the conclusion. 



2. Governing Equations 

We investigate the Parker instability with respect to how cosmic rays affect the fluid 
element in a uniformly-rotating galactic disk, under the influence of external gravity from 
the galactic center. A local rectangular domain representing the corotating sheet of the disk 
in the vertical plane is used for the simulations. The horizontal component of the radially 
inward external gravity is balanced by the centrifugal force, and hence, in the momentum 
equation, only the vertical component of external gravity and the Coriolis force are present 
to account for the rotation. Self-gravity from the plasma gas is not included. 

A hydrodynamic approach is adopted, in which the cosmic rays and the thermal plasma 
are two separate fluid components of a plasma system. The plasma fluid has a mass density 
of p, a thermal pressure of P g , and a cosmic-ray pressure of P c , all of which are threaded 
by a frozen-in magnetic field B. The cosmic-ray energy equation is described based on 
the diffusion-convection equation (Drury & Volk 1981; Jones & Kang 1990, Ko 1992), which 
treats the cosmic rays as a hot massless fluid and neglects the momentum spectrum of cosmic 
rays to simplify the governing equations. The artificial separation of the cosmic rays from 
the plasma helps to distinguish the role played by high- and low-energy components. 

The cosmic-ray diffusion-convection equation supplements the standard set of ideal 
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MHD equations. The governing equations in the complete vector form are written as 
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where V denotes the plasma fluid velocity; I denotes a unit tensor; 7 ff denotes the adiabatic 
index, i.e. ratio of the heat capacity at a constant pressure to that at a constant volume, of 
the thermal plasma gas; 7 C refers to the adiabatic index of the cosmic rays; k represents 
the cosmic-ray diffusivity; f2 represents the rotation angular frequency; and g is the external 
gravitational acceleration. Deriving this equation set involves use of the distribution function 
in Vlasov equation (Skilling 1975a, b, c). Such governing equations have also been used in 
the 2D hyperbolic tangent temperature model of Kuwabara, Nakamura & Ko (2004), and 
the 3D model of Lo, Ko & Wang (2011). 

In addition to balancing the energy equations, the term cosmic-ray pressure P c affects the 
momentum equation Eq. (|2]). However, particles of cosmic rays do not interact with plasma 
directly; they interact with plasma via the magnetic field. On a microscopic scale, resonant 
scattering of Alfven waves keeps the cosmic rays nearly isotropically distributed everywhere 
with respect to the thermal plasma background. A situation in which the cosmic-ray pressure 
possesses a gradient VP C influences the motion of the plasma gas. Thus, in formulating the 
equations, the interaction of cosmic-ray particles with a thermal plasma is represented by the 
cosmic-ray pressure P c and its gradient. The transport of the cosmic-ray pressure is described 
by a macroscopic diffusion coefficient, k , an energy-weighted mean diffusion tensor, defined 
as 



K±)bb + K±S; 



(6) 
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where b denotes a unit vector along the B direction. 



The limit of this model is mainly that one must assume a priori knowledge of the cosmic- 
ray pressure and energy density that satisfies the adiabatic index j c = 1 + P C /E C . Here, Eqs. 
([T])-(j5]) are solved numerically in a local reference frame in Cartesian coordinates, whose 
center lies at a galactocentric radius R Q and orbits the galaxy with a fixed angular velocity 
Q = Q(R ). In this local frame, the radial, azimuthal, and vertical spatial coordinates are 
related to the Cartesian coordinates such that x = R — R , y = R a ■ (4> — Q t), and z — z. In 
two dimensions a slice of flow is simulated, and so the rotational effect is not present. 



A finite difference method with operator splitting is employed to solve the governing 
equations. Eqs. (OJ)— Q are in the flux-conservative form and solved by 2-Steps Lax-Wendroff 
explicit method. The cosmic-ray energy equation Eq. ^ is divided into the convection and 
diffusion parts. The convection part is first converted to a conservative form and then 
also solved by 2-Steps Lax-Wendroff explicit method. The diffusion part is solved by the 
Bi-Conjugate Gradients Stabilized (BICGStab) implicit method. 

In the 2-step Lax-Wendroff method, Eqs. ([TJj-Q and the convection part of Eq. ^ 
are rewritten in the conservative form: 



where U can be density p, momentums pV, energy E g and magnetic field B; F,G,H are the 
flux of U in x-,y-,z- direction; and S is the source term. The x component of Eq. ([7]) is split 
into two steps: 



3. Numerical Algorithms 




(7) 
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where superscript n denotes time advection and subscript i represents space grid. 



The diffusion part of Eq. d5J) is written as: 




V ■ P?V(fi c )] = 0, 
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whose finite-difference form is expressed as Ax = b, where the submatrix in A = [A 1; A 2 , A 3 , 
A4,A 5 , A 6 , A 7 ] and the vector b is the function of time step At, the cosmic-ray diffusion 
coefficient k , the magnetic field B and the cosmic-ray energy at the nth iteration. The 
vector x is the unknown: 
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In order to solve x in Eq. (|9]), the BICGStab method is employed to handle asymmetric 
linear systems and reduce the operations per iteration to 0(N 2 ), where N is the number of 
unknowns in the discretized domain. This method is more efficient than the direct solution 
methods such as LU decomposition, which require 0(N 3 ) operations. The BICGStab is 
an iteration method that uses an initial guess x° to find a corresponding residual r° = 
b — Ax°, and then iterate to the i-th step r l to an accepted value by means of bi-conjugate 
matrix-vector operation. For comparison, when solving the nonlinear system F(y) = with 
Newton's method: 



Vi = Vi-i 



Vi-i + s h 



i-lj 



the correction 5i is determined by solving the gradient of the function, whereas in the 
BICGStab method, 



ro + 



\Vi-l\ 



is used for the iteration. The BICGStab treatment for cosmic-ray energy diffusion equation 
is an implicit method, and thus the Courant condition is not affected. 



4. Normalization and Initial Values 

The parameters and variables used in the governing equations withstand an extremely 
large contrast when expressed in dimensional units, possibly incurring significant errors due 
to the presence of ill-conditioned matrices and thus infeasible for numerical calculations. To 
overcome such numerical difficulties, the above equations are normalized to non-dimensional 
values that are close to unity. The non-dimensional values is denoted by the superscript 
'*' while the scaling factors are denoted by a '0' subscript. Some of the scaling factors also 
represent the quantities in equilibrium at the midplane of the Galactic disk. 
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The length variables are scaled based on the pressure scale height H Q : x* = x/H , 
y* = y/Ho, and z* = z/H . H is determined by integrating the hydrostatic equilibrium 
equation: 



d 
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given the ideal gas law, the values of sound speed C s o, the gravitational acceleration g z in 
the z-direction, the ratio of plasma pressure to magnetic pressure a, the ratio of plasma 
pressure to cosmic-ray pressure /3; and 7 S . The gas at the Galactic disk plane is nearly 
isothermal such that 7 3 ~ 1.0. The sound speed at the disk plane is C s q = 10 km/s. Thus, 
for a = (3 = 1 and g z = 2.28 x 10~ 8 cm s -2 , H = 50 pc = 1.54 x 10 20 cm. Given C s0 , the 
velocity is normalized to V* = V/C s0 , and the time scaling factor r is defined as H /C s0 , 
i.e. t = 1.54 x 10 14 sec ~ 5 Myr. The ISM density is used to normalize the density, with p 
being about one atom per cubic cm in the ISM, or po — 1-6 x 10~ 24 g cm -3 . 

With the above scaling parameters, the mass equation Eq. ([T]) becomes 
d(p p*)/d(r t*) + 1/HoV* ■ (p p*C s0 V*) = 0, 

which can be converted into dp* /<9t*+V*-(p*V*) = 0. Similarly, the factor p C s0 /(H /C s0 ) = 
PqCIq/ 'H at the left hand side of the momentum equation Eq. ^ is balanced by the same 
factor at the right hand side. The normalized plasma gas pressure is P* = P g /P g o, where 
Pgo — C*soPo — 1-6 x 10 -12 g cm -2 s~ 2 , and so the gradient of plasma pressure becomes 

Po c%v*p;/h . 



Scaling of the cosmic-ray pressure also employs the gas pressure, P* = P c /P g o- The 
scaling factor for the magnetic field, B , is determined by Bq = P g0 , and for the Galaxy, 
B = ^/P^ = 1.26 x 10~ 6 Gauss. 



The gravitational acceleration is scaled to go = C 2 /ifo — 6.49 x 10 -9 cm s -2 , while 
the scaling factor for the rotating angular frequency is the reverse of time: Qq — C s q/Hq = 
0.65 x 10 -26 s _1 . Using the same scaling parameters, the momentum equation Eq. ([2]), the 
induction equation Eq. ^ and energy equation Eq. Q can all be converted to dimensionless 
values. 

Finally, in Eq. ([5]), the cosmic-ray diffusion coefficient k is normalized to k* = k/kq, 
where Kq = C s qHq = 1.54 x 10 24 cm 2 s _1 . Theoretical calculations using path length distri- 
bution, life-time of radioactive secondary CR, etc. suggested that the diffusion coefficient for 
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CR above 1 GeV in the Milky Way Galaxy is 1 - 3 x 10 28 cm 2 s" 1 (Berezinskii et al. 1990; 
Ptuskin 2001). Due to numerical difficulties with strong magnetic field and large diffusion 
coefficient, this work adopts a magnetic field strength of several pG, like those in the solar 
vicinity, and a ~ 10 — 100 times lower k than suggested. Crocker et al. (2010) suggest a 
lower limit of 50 pG on 400 pc scales near the Galactic center. However, measurements of 
the amplitude of the magnetic field depend on the spatial scale and the energy equipartition 
or pressure equilibrium among various Galactic components, which may differ by 2 orders of 
magnitude. Pure MHD models with strong magnetic field have been presented in Machida 
et al. (2009) and Takahashi et al. (2009). 



5. Initial Equilibrium Background 

We adopt two initial equilibrium backgrounds: a temperature distribution that follows 
a hyperbolic tangent profile and a density distribution that follows an exponential profile. 
Both these backgrounds bear a declining density profile with height. They differ mainly 
in that the hyperbolic tangent temperature has an ascending transition zone that divides 
the distribution into two distinct regions, while the exponential density follows a continuous 
decline. Given the temperature or density profile, the profiles for the density /temperature, 
pressure, magnetic field and other variables can be derived from the hydrostatic condition. 

In the isothermal case the criteria for Parker instability to develop is -^(B/p) < 0. 
When the rising field lines grow, the flow becomes unstable. In our cases, the growth of the 
instability depends on additional parameters such as the width of transition region and the 
height of Galactic halo (disk thickness). 



5.1. Hyperbolic Tangent Temperature 



The hyperbolic tangent temperature profile is a two-temperature, layered disk (Shibata 
et al. 1989a), described by 



C*(z) = T(z) = T + (T h -T 



tanh 



W tr 



+ 1 



(12) 



where C s (z) is the sound speed, T = 10 4 K is the disk temperature and T h = 25T is the 
halo temperature. The initial dimensionless temperature is 1.0, equivalent to 10 4 K. Given 
the ideal gas law P g = pT/jg and the gravitational acceleration g z , the initial density profile 



is solved using the hydrostatic equation Eq. (10). A constant acceleration g z is assumed 



because the grid domain is small and not far way from the disk plane. The dimensionless 
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gravitational acceleration is set as 1-0 /j g (dimensionless units), equivalent to 6.49/7 9 x 
10~ 9 cm s~ 2 , or ~ 6.18 x 10~ 9 cm s~ 2 , for ^ g = 1.05. Acceleration can vary widely from 
galaxy to galaxy. Although a little higher, this value is comparable to the value derived from 
observations of the spatial density distribution and the velocity-distance correlation 500 pc 
above the Galactic disk midplane, 4 x 10~ 9 cm s -2 (Oort 1960; Bahcall 1984; Kuijken & 
Gilmore 1989). The initial density at the disk plane is 1.0 (dimensionless units), equivalent 
to 1.6 x 10~ 24 g cm -3 . 

After obtaining the density profile, the plasma pressure profile P g = C^p/-y g is derived. 
The magnetic field is assumed to align in the x-direction and vary with height z, and B x (z) = 
^jRitPg/a. The initial dimensionless pressure at the disk is 1.0 (2 x 10~ 12 g cm -2 s~ 2 ). The 
initial magnetic field at the disk for a = 1 is \/8n ~ 5 (dimensionless units), or 5B = 6.34 
(//Gauss), This field strength is very close to the radio synchrotron measurement of 6 //Gauss 
averaged over a radius of about 1 kpc around the Sun (Beck 2009). 

The initial cosmic-ray pressure at disk for = 1 is 1.0 (dimensionless units), or 2 x 
10~ 12 g cm" 1 s~ 2 . The rotating angular frequency Q (not used in this work) is a free 
parameter; Q* = 1 gives 0.65 x 10~ 26 s -1 , or approximately 7 times the angular frequency 
at our Sun, Q & = 220/7.6 km s _1 kpc -1 . The initial equilibrium background is assumed 
quiescent. No initial y- velocity and no rotation effects are present in the system. 



5.2. Exponential density 

The exponential density case produces a nearly isothermal background that is similar 
to the isothermal model of Kim et al. (1998). However, since 7 9 = 1.0 imposes singularity 
in the energy equation, 7 g > 1 is adopted. The exponential density profile is defined as: 

p = p exp y-jfj ' ( 13 ) 

where po = 1.0 (dimensionless units), equivalent to 1.6 x 10~ 24 g cm -3 ; the density scale 
Hp = 4.0 (dimensionless units) is equivalent to 6.16 x 10 20 cm. The plasma pressure follows 



P g = pT/^f g . Using Eq. (10), the temperature becomes: 



Tt= <pfcnl' (14) 

which is a constant. A minimum value of the adiabatic index 7 9 = 1.01 is adopted in this 
model. Similar to the hyperbolic tangent model, the initial magnetic field and cosmic-ray 
pressure are defined through the parameter a and f3, B x (z) = y^8nPg/a and P c = P g //3. 
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6. Grid Setup and Perturbation 



2D simulations are performed in a rectangular domain in the x — z plane. To excite the 
instability, velocity perturbations are added to the quiescent background. Two perturbation 
forms are examined: eigenmode sinusoidal wave and random seed. The strong explosive 
perturbation in CR pressure developed from supernova explosions has been considered in 
Kuwabara et al. (2004). Notably, although supernova shocks has been believed as efficient 
cosmic ray accelerator, channeling at least 10% of the supernova explosion energy into cos- 
mic rays, signature of cosmic-ray protons, the key observational evidence in support of the 
claim that isolated supernova remnants (SNRs) are the main accelerators of galactic CRs, 
has remained elusive after decades of observational effort (Allen, Houck, & Sturner 2009). 
Calculations of the maximum energies of shock-accelerated electrons (Reynolds & Keohane 
1999) and hydrodynamic simulations of Rayleigh- Taylor instabilities in SNRs (Ferrand et al. 
2010, Fraschetti et al. 2010, Wang 2011) both indicate that young SNRs could not be re- 
sponsible for the highest-energy Galactic cosmic rays, unless an unrealistically high injection 
rate of cosmic rays, and thus an enormous energy loss in SNRs, is invoked. Precise measure- 
ments of cosmic-ray proton and helium spectra by Adriani et al. (2011) furher excludes a 
single power law energy distribution for the two species and thus challenged the supernova 
acceleration paradigm. 

The sinusoidal perturbation employs a perturbing velocity V x described by 



where A = 20H is the most unstable wavelength derived from the linear analysis for k» = 
200 (Kuwabara et al. 2004). V x is applied within the region of 4Hq < z < 8Hq and 
\x — xo\ < 1/2A, where xq = 40 (dimensionless unit). The grid consists of 102 x 301 
zones in a rectangular region of 4 kpc x 3.5 kpc. The horizontal length of each grid zone is 
Ax = 0.8H , while the vertical length is nonuniform: 



Az is increased above z = 25H by a ratio of 1.05, until Az max = 5 A z (at z = 0). The 
upper grid domain accommodates the unperturbed magnetic field to avoid outflows across 
the grid boundary and spurious results. 

For the random velocity perturbation (Shore & Larosa 1999), a series of random veloc- 
ities V x and V z with a maximum amplitude of 5% is added horizontally into the background 
(Chou et al. 2000). The grid covers an area of 9 kpc x 11.5 kpc using 512 x 512 zones, with 





0.15# , < z < 25#< 

min{ 1. 05 Az, Az max }, otherwise 



o 
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Ax = 0.35H and 

Az _( 0.15-Ho, 0<z<35H 

\ min{1.05Az, Az max }, otherwise 



7. Results 

7.1. Pure MHD Models 

We examine how cosmic-ray diffusion affects Parker instability by comparing a series 
of pure MHD simulations with simulations that include the cosmic-ray effect. Figure 1 
illustrates those results using sinusoidal eigenmode perturbations. Figure 1(a) presents the 
hyperbolic tangent temperature model, while Figure 1(b) displays the exponential density 
model at the same epoch. In these figures a grid domain of z < 60 that just lies within 
the maximum height of the magnetic loops is shown. Logarithmic values of variables in a 
dimensionless unit are presented in color contours. White solid lines represent the magnetic 
field lines, and small white arrows refer to the velocity vectors. The superscript '*' denoting 
dimensionless values is omitted hereafter. The exponential model shows a faster growth of 
the instability because the decline in the gas density facilitates the instability. 



7.2. Cosmic-Ray Diffusion Coefficient 

Figure 2 compares the maximum height of the magnetic loops versus the cosmic-ray 
diffusion coefficient k\\ in both models at t — 33, when the numerical time step becomes very 
restrained. A smaller CR diffusion coefficient means stronger coupling between the cosmic 
rays and the plasma. A maximum of k\\ = 200 (1.54 x 10 24 cm 2 s" 1 ) is examined. In 
both models, the instability becomes more prominent with an increasing k\\ , and the mixing 
extends to a similarly maximum height of z ~ 55. Because the gradient of cosmic-ray pressure 
VP C declines with an increasing diffusivity, cosmic-ray diffusion facilitates flow instability. 
Kuwabara et al. (2004) studied the instability in the hyperbolic tangent temperature case. 
For small-Kii, the instability is impeded by the CR pressure gradient force that interferes 
with the falling motion of the matter, while for large-Ky, the magnetic loop can rise to larger 
scales. Our simulations reach a similar outcome. In Kuwabara et al (2004), the case of 
strong explosive perturbation in CR pressure presumably from supernova yielded a larger 
growth rate for a smaller diffusion coefficient, but such a reversed correlation only occurs in 
the early stage; in the later stage the growth rate becomes smaller when compared to that 
of a large diffusion coefficient model. In the 3D models of Lo, Ko & Wang (2011) employing 
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sinusoidal perturbations, however, the instability growth first rises with k\\ but then drops 
when k\\ > 32, because under strong diffusion, the falling of the gas becomes supersonic; as 
a result of shock, the cosmic rays are redistributed with a smoother pressure gradient, which 
tends to stabilize the flow. 



7.3. Parameter a and [5 

The scale height H = 50 pc is set for a — f3 — 1. Varying a = P 9 /Pb modifies the length 
scale, resulting in different flow structures. Figure 3 presents snapshots of the hyperbolic 
tangent temperature model with a tenfold gas pressure, a = 10.0; the flow extends to z ~ 40, 
which displays more prominent mixing than the a — 1 case [z ~ 20). Other parameters used 
are (3 = 1.0; 7 S = 1.05; 7 C = 4/3; Wt r = 0.6 (dimensionless units, 30 pc or 0.92 x 10 20 cm; 
z h = 18 (900 pc or 1.39 x 10 23 cm; rc„ = 2.0 (1.54 x 10 23 cm 2 s" 1 ); and k ± = 0.0. 

7.4. Parameter 7 S and w tr 

The temperature distribution is assumed to be uniform in the exponential density model. 
In this case, modifying a and (5 changes only the isothermal temperature value. Notably, 
the parameter (3 = P g /P c should also influence the scale length (see the scaling for H ). 
However, because cosmic rays largely diffuse along the magnetic field lines, varying (3 does 
not significantly affect the scale height. 

The isothermal case is approximated by using a minimum adiabatic index of 7 S = 1.01. 
In an ISM, ^ g should be significantly smaller than the ideal value 5/3 and close to 1.0, 
because thermal instability smooths out the temperature gradient. However, because the 
diffusion of cosmic rays along the magnetic field lines is only slight and the instability growth 
expedites for more condensed gas, thermal instability is expected to be deterred when cosmic- 
ray diffusion is present (Shadmehri 2009). This phenomenon may resemble the effect of 
increasing 7 9 above the typical value 1.05 that is used in most of the simulations for Parker 
instability. According to our results, instability is dampened as 7 increases; with 7 = 5/3, 
the flow remains quiescent at t ~ 98 (Figure 4a). 

In the hyperbolic tangent temperature model, an increasing width of transition region 
w tr makes the instability less feasible (Figure 4b). Suppression of the instability is attributed 
to the rapid rise of temperature and the flattening of density near the transition region. 
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7.5. Random Perturbation 

The case using random velocity perturbation may reflect the circumstances closer to the 
actual situation in an ISM. In this case, the development of instability is slower than the case 
using eigenmode perturbation, Figure 5 reveals that the exponential density model exhibits 
slightly higher floatation of the magnetic loops, but similarly strong mixing and filamentary 
structures are induced in both models. 



7.6. Trebly Sinusoidal Perturbation with kj_ 

The perpendicular or cross field lines diffusion coefficient k± (see Eq. [6]) is often sub- 
stantially smaller than the parallel coefficient Ki\, which is only around 2% — 4% of ku (Ryu 
et al. 2003). We found that when incorporating cosmic-ray cross field line diffusion k±, the 
magnetic field lines are aligned more vertically to the disk. As such a process may be re- 
lated to the acceleration of galactic wind, a new initial condition using the same equilibrium 
background is invoked to examine the case involving kj_- 

A grid of 512 x 512 zones in a domain of 9kpc x 11.5kpc is adopted, with each zone 
having a size identical to the case of random perturbation. The applied velocity perturbation 
is similar to our previous sinusoidal eigen mode perturbation, but within a finite rectangular 
region of \x — Xq\ < 3/2A, where Xq = 90.0, along with other similar parameters: w tr = 0.6; 
z h = 18; a = $ = 1; j g = 1.05; 7 C = 4/3; k\\ = 2.0; and k± = 0.02. 

Figure 6 shows the exponential case at the dimensionless time t ~ 30 (150 Myr). Three 
equally-sized loops arise and extend to z > 50. The flow pattern resembles the hyperbolic 
tangent case in Kuwabara et al. (2004). The interaction between loops becomes more 
pronounced at later times; the central one becomes compressed by two outer adjacent loops, 
while the loops continually rise to a height of 2.0 kpc. At t — 60, the flow reaches a maximum 
vertical height and becomes saturated. Different from the case without the flow velocity 
is significantly increased such that the loops rise to a much higher altitude, z > 80. Therefore, 
even at a value of 1% of K\\, the effect of k± on the mixing appears significant (Figure 6). 
Although such a result contradicts the estimate made by linear analysis (Ryu et al. 2003), 
it indicates that Parker instability is likely to contribute to the galactic wind acceleration 
perpendicular to the galactic disk (Wiegelmann, Schindler, & Neukirch 1997). 
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8. Conclusion 

This study has elucidated the development of Parker instability undergoing cosmic-ray 
diffusion using 2D MHD simulations. Exactly how the initial conditions affect the growth of 
instability is examined. Two equilibrium backgrounds are constructed based on hyperbolic 
tangent temperature and exponential density, with variations in the physical parameters and 
perturbations. 

Our simulations gives rise to a structure that is closely resembles the long, thin wave-like 
sheared helmet plumes in the solar corona, caused by magnetic buoyancy (Shore & Larosa 
1999). The perpendicular component of cosmic-ray diffusion coefficient k± is included as 
well. When k± is only l% — 2% of the parallel component rey, despite the small development, 
vertical magnetic structures and outward flow arise, particularly in the exponential density 
model. The substantial increases in the velocity and altitude of the unstable flow reveal 
the importance of cosmic rays to the galactic wind. In general, in the exponential density 
model, the morphology of the clumpy structure is filamentary. Conversely, in the hyperbolic 
tangent temperature model, a more concentrated and round shape of clumps like the giant 
molecular cloud are observed at the foot points of rising magnetic arches. The growth of 
instability in the hyperbolic tangent model is less rapid than that of the exponential density 
case since the pressure in the exponential density decreases faster with an increasing height. 

The small gradient of cosmic ray pressure along the z direction caused by diffusion 
explains why an increase in k\\ facilitates the growth of the unstable undular mode. As a 
result, the gas falls down more rapidly, and adjacent loops join together to form a large 
loop-like bubble. 

Exponential density with the adiabatic index 7 S ~ 1 (i.e. isothermal) produces the 
highest magnetic loops and most flow instability. A decreasing 7 S leads to more condensed 
gas and, ultimately, a more unstable flow. While examining a situation without cosmic- 
ray diffusion, Parker indicated that the criterion for instability is 7 9 < 1.36. Although the 
undular mode is expected to be suppressed when ^ g > 1.4, numerical results in our cases using 
1.3 < 7 9 < 1.4 indicate that the high adiabatic index still produces a non-uniform density 
distribution before the simulations end due to stringent numerical conditions. Therefore, the 
instability may still occur for large 7 if a stronger perturbation such as a supernova explosion 
is invoked. 

The hyperbolic tangent temperature model yields a density distribution that increases 
with height, whose steepness depends on the width of the transition zone w tr . For a smaller 
w tr , instability is increased as the density gradient correspondingly diminishes near the 
unstable regions. 
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When a larger ratio of plasma pressure to magnetic pressure a is applied, the pressure 
inside the flux tube diminishes, thus the instability growth is facilitated. Also, as the length 
scale decreases accordingly, unstable structures are observed on a smaller scale. Contrarily, 
a small a yields a larger scale height and so complicates the mixing of an unstable flow. 

Finally, the ratio of plasma pressure to cosmic-ray pressure j3 does not influence the 
length scale since the diffusion of the cosmic-ray pressure does not contribute to the scale 
height. 
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Fig. 1. — Logarithmic density contours in the pure MHD case developed from a sinusoidal 
perturbation with a — 1 at t ~ 69. (a) A case of hyperbolic tangent temperature; (b) a 
case of exponential density. 
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Fig. 2. — Maximum height of loops z versus cosmic-ray diffusion coefficient k\\ at t ~ 33. 
Diamond denotes the exponential density model and square denote the hyperbolic tangent 
temperature model. 
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Fig. 3. — Snapshots of logarithmic values of density, pressure, and cosmic-ray pressure 
contours for hyperbolic tangent temperature model at t ~ 33 with a = 10 and k± = 0. 
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Fig. 4. — Density contours with k± = for (a) 7 9 = 5/3 and (b) w tr = 10.0 at t ~ 98 in the 
hyperbolic tangent model. Figure (a) shows small fluctuations while Figure (b) does not. 
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Fig. 5. — Density contours with k± = at t ~ 49 for (a) hyperbolic tangent temperature 
(b) exponential density, using random perturbation. 




Fig. 6. — Density contour of the exponential density model at t ~ 32 using trebly sinusoidal 
perturbation and a cosmic-ray cross field line diffusion k,j_ = 2%. The high velocity indicates 
Parker instability may contribute to the galactic winds. 



